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Abstract. The dynamics, time evolution of the mass dis- 
tribution, and gravitational wave signature of coalesc- 
ing neutron stars described by polytropes are compared 
with three simulations published previously: (a) "Run 2" 
of Zhuge et al. (1994), (b) "Model III" of Shibata et 
al. (1992), and (c) "Model A64" of Ruffert et al. (1996). 
We aim at studying the differences due to the use of differ- 
ent numerical methods, different implementations of the 
gravitational wave backreaction, and different equations 
of state. We integrate the three-dimensional Newtonian 
equations of hydrodynamics by the Riemann-solver based 
"Piecewise Parabolic Method" on an equidistant Carte- 
sian grid. 

Comparison (a) confronts the results of our grid-based 
PPM scheme with those from an SPH code. We find that 
due to the lower numerical viscosity of the PPM code, the 
post-merging oscillations and pulsations can be followed 
for a longer time and lead to larger secondary and ter- 
tiary maxima of the gravitational wave luminosity and to 
a stronger peak of the gravitational wave spectrum at a 
frequency of about / » 1,8 KHz when compared to the 
results of Zhuge et al. (1994). In case (b) two grid based 
codes with the same backreaction formalism but differ- 
ing hydrodynamic integrators and slightly different initial 
conditions are compared. Instead of rotationally deformed 
initial neutron stars we use spherically shaped stars. Sat- 
isfactory agreement of the amplitude of the gravitational 
wave luminosity is established, although due to the differ- 
ent initial conditions a small time delay develops in the 
onset of the dynamical instability setting in when the two 
stars come very close. In (c) we find that using a polytropic 
equation of state instead of the high-density equation of 
state of Lattimer & Swesty (1991) employed by Ruffert 
et al. (1996) does not change the overall dynamical evo- 
lution of the merger and yields agreement of the gravi- 
tational wave signature to within 20% accuracy. Whereas 



the polytropic law describes the dynamical behaviour of 
the bulk of the matter at and above nuclear density suf- 
ficiently well, we, however, find clear differences of the 
structure and evolution of the outer layers of the neutron 
stars where the stiffness of the equation of state is largely 
overestimated. This has important implications for ques- 
tions like mass loss and disk formation during the merging 
of binary neutron stars. 
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1. Introduction 

Coalescing neutron stars are interesting and powerful 
sources of gravitational waves. These mergers — together 
with black holes, rapidly rotating neutron stars, and su- 
pernovae — are among the most promising candidates 
(Thorne 1992) to be detected with the large gravita- 
tional wave experiments that will soon be operational or 
are currently planned, LIGO (Abramovici et al. 1992), 
VIRGO (Bradaschia et al. 1991), GEO600 (Danzmann 
et al. 1995a), and the space-based LISA (Danzmann et 
al. 1995b). The detectability does not only depend on the 
strength of the emission, but also on the frequency of the 
gravitational waves relative to the resonance frequencies of 
the detectors (Finn & Chernoff 1993; Finn 1994). Only an 
a priori knowledge of the structure of typical gravitational 
wave signals allows the extraction of faint signals from a 
very noisy background. Theoretically predicted wave tem- 
plates have to be produced by numerical simulations. Ad- 
ditionally to being sources for gravitational waves, coalesc- 
ing neutron stars (and its relative, the merger of a neutron 
star with a black hole) remain the favorite cosmological 
model for gamma-ray bursters (Dermer & Weiler 1995) 
because the sources are known to exist, the right amount 
of energy is involved, and the expected merger rates are 
in agreement with the required gamma-ray burst rates. 



2 



M. Ruffert et al.: Coalescing neutron stars - gravitational waves from polytropic models 



The search for gravitational waves thus calls for mod- 
els of coalescing neutron stars that are numerically highly 
resolved and as realistic as possible. Several groups have 
attempted such investigations and have published results. 
In a series of papers starting with Oohara & Nakamura 
(1989) a Japanese group used variations of a grid-based 
code to find the gravitational wave luminosity for differ- 
ent conditions, e.g. different masses, initial separation or 
spins of the neutron stars (Shibata et al. 1993, and ref- 
erences therein). Zhuge et al. (1994) concentrated on the 
energy spectrum that gives the energy emitted in gravita- 
tional waves at different frequencies. They discussed the 
information about the dynamics of the coalescence that 
may be extracted from a measured spectrum. Focussing 
on dynamics and stability, Rasio & Shapiro (1994) found 
that triaxial configurations are formed by the merging of 
neutron stars with a sufficiently stiff polytropic equation 
of state. Non-axisymmetric structures occur for adiabatic 
exponents of 3 or higher (polytropic index 0.5 or below). 
In this case the peak amplitude of the gravitational wave 
emission is substantially larger and the emission proceeds 
for a longer time after the coalescence. 

Our treatment of the gravitational wave backreaction 
is identical to what has been used by Shibata et al. (1992). 
We consider this procedure to be more accurate than the 
point-mass quadrupole approximation used by Davies et 
al. (1994) and Zhuge et al. (1994). Their approximation 
was used when the neutron stars were still separated and 
subsequently switched off when the merging took place. 
Although Rasio & Shapiro (1994) investigated the stabil- 
ity of the neutron star binary orbit in detail, they did not 
take into account the gravitational radiation reaction at 
all. 

The three investigations done by Davies et al. (1994), 
Zhuge et al. (1994) and Rasio & Shapiro (1994) all used 
SPH codes, whereas our simulations, as well as the ones 
by Shibata et al. (1992), were performed with explicit, 
Eulcrian, finite-difference grid-based schemes. The algo- 
rithm we use is based on the PPM-scheme of Colella & 
Woodward (1984) which is a higher-order Godunov-type 
method in that it locally solves Riemann problems. It is 
thus superior to a particle-based method like SPH in the 
handling of shocks and also has a rather small numerical 
viscosity compared to other Eulerian algorithms and SPH. 

We chose to compare our numerical scheme with the 
following three calculations. (1) "Run 2" of Zhuge et 
al. (1994) was chosen since in this model the neutron 
stars have masses and radii similar to those in Ruffert 
et al. (1996) (although the equations of state were differ- 
ent). Additionally, the published detailed spectral analysis 
of the gravitational wave emission is well suited for com- 
parisons. (2) "Model III" of Shibata et al. (1992) was se- 
lected, because their numerical scheme is grid-based and 
the algorithm to treat the gravitational wave backreaction 
is similar to ours. Therefore differences of the results are 
associated with differing hydrodynamical integrators and 



numerical resolution. In both cases (1) and (2), polytopic 
equations of state were used which we also applied in our 
comparative calculations. (3) "Model A64" of Ruffert et 
al. (1996) served for comparison in order to investigate the 
differences caused by the use of a polytropic equation of 
state instead of the more realistic equation of state of Lat- 
timer k Swesty (1991) which allowed Ruffert et al. (1996) 
also to include the effects due to neutrino emission. 

We did not select a model computed by Davies et 
al. (1994) because although their computations were done 
with a polytropic equation of state, they did not give infor- 
mation about the gravitational wave emission, since their 
formula was too crude. Rasio & Shapiro (1994), on the 
other hand, were very careful in constructing equilibrium 
models and in investigating the stability of these models 
as well as the gravitational wave forms and luminosity. 
Nevertheless, we preferred for our present investigations a 
comparison with a run of Zhuge et al. (1994), because (a) 
Rasio & Shapiro (1994) did not include any backreaction 
of the gravitational waves, and (b) the spectra presented 
by Zhuge et al. (1994) were better suited for comparisons. 

Section gives a brief description of the numerical 
methods, introduces the set of computed models and de- 
tails the initial conditions that our simulations are started 
from. The results of our computations and the compari- 
son with previous published models are presented in the 
three Sects. ||, |], and ||, separately. A summary follows in 
Sect. | 

2. Computational procedure and initial conditions 

The numerical procedures implemented to simulate our 
models of coalescing neutron stars were described in detail 
in a first paper (Ruffert et al. 1996). In the following, we 
shall specify only the differences of the treatment of the 
presented polytropic models compared to what has been 
described in Ruffert et al. (1996). 

2.1. Hydrodynamics and equation of state 

The neutron stars are evolved hydrodynamically via the 
Piecewise Parabolic Method (PPM) developed by Colella 
and Woodward (1984). The code is basically Newtonian, 
but contains the terms necessary to describe gravitational 
waves and their backreaction on the hydrodynamical flow 
(Blanchet, Damour & Schafer 1990) via energy and an- 
gular momentum loss. The volume integral over the local 
emission and backreaction terms yields the total gravita- 
tional wave luminosity. All spatial derivatives necessary to 
compute the gravitational wave terms are implemented as 
standard centered differences on the grid. 

The Poisson equations in integral form, necessary for 
the calculation of the gravitational potential as well as for 
the wave backreaction, are interpreted as convolution and 
calculated by fast Fourier transform routines: non-periodic 
grid boundaries are enforced by zero-padding (e.g. Press 
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Table 1. Parameters of the three Models ZCM, SNO, and RJS. R is the initial radius of the two identical neutron stars, po their 
initial central density, ao the initial center-to-center separation of the stars, M the mass of the neutron stars, K the polytropic 
constant, F the adiabatic exponent, L the grid size, N the number of zones per dimension, if the total time interval of the run, 
£ the total energy emitted in gravitational waves within time if, C is the maximum gravitational- wave luminosity, and rh is the 
maximum amplitude of the gravitational waveforms observed from a distance r 
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0.88 


Zhuge et al. (1994), Run 2 
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13.6 


Shibata et al. (1992), Model III 
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15 


5.7 


42 


1.63 


2.438 


2.319 


82 


64 


4.6 


2.8 


2.6 


2.07 


Ruffert et al. (1996), Model A64 



et al. 1992; Eastwood & Brownrigg 1979). The accelera- 
tions that follow from the potential are added as source 
terms into the PPM algorithm. For models ZCM and RJS 
the number of cubic grid zones in the orbital plane is 
64 by 64, while perpendicular to the orbital plane 16 zones 
were used with the same zone size in all three dimensions. 
We use double the amount of zones (128) in every spatial 
direction for model SNO in order to have a similar number 
of zones (121) as Shibata et al. (1992). We take advantage 
of the symmetry about the orbital plane to actually cal- 
culate only one side of the plane. 
A polytropic equation of state, 

P = Kp T , (1) 

is used which relates the density p to the pressure P via 
two parameters, the polytropic constant K and the adi- 
abatic exponent T = 1 + i, with n being the polytropic 
index. This equation with K and T chosen to be the same 
in the whole star was employed to construct the initial 
model of the neutron star. Since the PPM code is writ- 
ten in conservative form (not taking into account source 
terms like gravity), the primary variables are the density, 
the momentum, and the total energy. One therefore has 
the choice of using either Eq. [I] to calculate the pressure 
from the density or the relation 

P=(r-l)e (2) 

to deduce the pressure from the internal energy density 
e. Analytically these formulations are equivalent, but nu- 
merical errors may destroy this equivalence. When shocks 
or viscous shear dissipation occurs in the hydrodynamic 
calculations, the use of Eq. || is preferable because due to 
the dissipative processes the parameter K of Eq. ^ is not a 
strict constant in space and time during the simulations. 
If one postulates analytical equivalence of the pressure 
in Eq. || with the polytropic equation of state, Eq. |], the 
generation of thermal energy in shocks manifests itself in 
a variation of K. Eq. [j] was used to construct the initial 
model and Eq. ^ to compute the hydrodynamic evolution. 
Due to these simple equations of state, effects of neutrino 



emission, absorption or annihilation were not included in 
the present simulations. 

2.2. Initial setup 

Our initial setup involves placing two identical spherical 
neutron stars at some distance from each other and giv- 
ing them an initial velocity. We prescribe the initial orbital 
velocities of the coalescing neutron stars according to the 
motions of point masses, as computed from the quadrupole 
formula. The tangential components of the velocities of the 
neutron star centers correspond to Kepler velocities on cir- 
cular orbits, while the radial velocity components reflect 
the emission of gravitational waves leading to the inspiral 
of the orbiting bodies (for details, see Ruffert et al. 1996). 
A spin of the neutron stars around their respective cen- 
ters is not added. Thus the initial state of the neutron star 
binary is defined by the following parameters: (a) center- 
to-center distance between the neutron stars and (b) three 
out of the following five quantities: mass, radius, and cen- 
tral density of the neutron star, and polytropic constant 
K, and adiabatic exponent T of the equation of state. 

The initial distance at which the neutron stars are 
placed at the beginning of our simulations is a compromise 
between physical accuracy and computational resources. 
On the one hand we would like to start with the neu- 
tron stars at a large distance in order to correctly sim- 
ulate the tidal deformation they experience during inspi- 
ral. However, the time to coalescence, which is propor- 
tional to the fourth power of the distance, increases by a 
huge amount if the distance between the neutron stars is 
only slightly increased. On the other hand, tidal deforma- 
tion studies (e.g. Lai et al. 1994; Reisenegger & Goldreich 
1994) show that tidal deformations for extended objects 
are only around 20% (for the principal axis) and normal 
mode excitation of at most a few percent of the neutron 
star radius at a distance of about 2.8 radii. 

For numerical reasons it is not possible to do the sim- 
ulations with arbitrary low density in the neutron star's 
environment. So the density of the matter in the surround- 
ings of the neutron stars is set to 10 9 g/cm 3 . To avoid 
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Fig. 1. Cuts in the orbital plane of Model ZCM at four instants in time showing the density contours together with the velocity 
field. The density contours are logarithmically spaced with intervals of 0.5 dex. The density is measured in units of g/cm 3 . The 
bold contours are labeled with their respective values. The legend at the top right corner of each panel gives the scale of the 
velocity vectors and the time elapsed since the beginning of the simulation 



this low-density matter being accelerated and falling onto 
the neutron stars, the velocities and kinetic energy are 
reduced to zero in zones in which the density is less than 
about 3-10 9 g/cm 3 . With a numerical grid of roughly 1 km 
spatial resolution one cannot accurately represent the den- 
sity decline near the surface of the neutron stars. We ar- 
tificially soften the edge by imposing a maximal density 
change of 2 orders of magnitude from zone to zone. The 
thickness of the surface layer thus results to about 3 grid 
zones. To avoid short numerical time steps, we changed the 
constant K of the equation of state in the 3 surface zones 
such that the sound speed was constant within this broad- 



ened surface layer. The same manipulation was applied for 
the surrounding gas at densities of about 10 9 g/cm 3 . 

We simulate three models with polytropic equations 
of state, the parameters of which are listed in Table El 
The initial conditions of each model are discussed in the 
respective section of the corresponding model. The cal- 
culations were performed on a Cray J90 8/512, needed 
about 3.2 MWords of main memory and approximately 
25 CPU-hours for models ZCM and RJS, and 27 MWords 
and 140 CPU-hours for model SNO. 
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3. Results of Model ZCM 



3.1. Initial conditions 

We chose Run 2 of Zhuge et al. (1994), because the mass 
and radius of the neutron stars are closest to those of the 
stars which we have simulated in Ruffert et al. (1996). For 
Run 2, moreover, Zhuge et al. (1994) show the gravita- 
tional wave forms, the gravitational wave luminosity, and 
the energy spectrum. This case yields a direct comparison 
between the SPH computations with ad hoc gravitational 
backreaction done by Zhuge et al. (1994) and our Eulerian 
calculation of Model ZCM where backreaction is included. 

Model ZCM has the same neutron star mass and ra- 
dius as well as an adiabatic exponent of T = 2 as used by 
Zhuge et al. (1994) for their Run 2. In order to save CPU 
time we place the neutron stars at an initial distance of 
45 km and not at the 60 km chosen by Zhuge et al. (1994). 
They were able to follow the inspiral from such a large dis- 
tance because their SPH simulation was done with a com- 
paratively small number of particles (N = 1024) which is 
computationally much cheaper than our grid-based sim- 
ulations. Unfortunately, they do not show a plot of the 
neutron star separation as a function of time for Run 2. 
Without this relation the connection between the tempo- 
ral evolution of their model and ours is not determined 
precisely. Thus we link our time to theirs by matching 
both calculations at the position of the first maximum of 
the gravitational wave luminosity (see Fig. ||) independent 
of how long it takes to spiral in from the initial separa- 
tion of 60 km (used by Zhuge et al. 1994) to our initial 
separation of 45 km. 

However, the information given by Zhuge et al. (1994) 
for their Run 1 can be used to interpret their Run 2. From 
their Fig. 10 one can infer at which time the separation of 
their neutron stars is 3i?, and then inspect their Fig. la to 
get an impression of the tidal deformations at that time: 
tidal bulges are visible but not dominant, in agreement 
with the analytical estimate of 20% deformation of the 



energy, model ZCM 



principal axis at a distance of 2.8R (sec Sect. 2.2). We did 
not include these bulges in our initial setup of the neutron 
stars. 

3.2. Dynamical evolution 

Snapshots of the density distribution in the orbital plane 
together with the velocity field are shown for four times in 
Fig. 0. These plots give an impression of how the merging 
process proceeds. Within less than one orbital period the 
merging is well in progress (Fig. 0a). The typical tran- 
sient spiral-arm pattern develops (Fig. 0b) due to tidal 
and centrifugal forces. One notices that for a fairly long 
time (for over 2 ms after t — 4 ms; e.g. Fig. 0c) the merged 
object has sort of a "double core" structure, visible by the 
dumb-bell shaped central contours. The quadrupolar de- 
formation seems to be less pronounced in the model of 
Zhuge et al. (1994). It is damped only slowly and indica- 
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Fig. 3. The gravitational wave luminosity as a function of time 
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lines 
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Fig. 4. The gravitational wave forms, h+ and h x , for 
Model ZCM 

tions of it are still present when we stop our simulation at 
t = 8.7 ms (Fig. gd). 

The temporal evolution of different energies is shown 
in Fig. |2| One notices that the coalescence is very smooth 
and that the potential energy decreases only by about 
15%. The internal energy first decreases a bit due to 
tidal stretching, then increases again because some ki- 
netic energy is dissipated into thermal energy, indicating 
the presence of shocks and friction. During the spiral-in 
and merging the maximum density decreases monotoni- 
cally from roughly 6-10 14 g/cm 3 to 5-10 14 g/cm 3 , because 
of tidal stretching before coalescence and internal "heat- 
ing" due to shear effects which dissipate kinetic energy 
during coalescence.El The total energy is conserved to 2% 
compared to the maximum potential energy. 

3. 3. Gravitational wave forms and luminosity 

The temporal structure of the gravitational wave luminos- 
ity is shown in Fig. ^ (the time integral of which can be 
seen in Fig. |^) . The luminosity during the main merging 
event is of similar size in our Model ZCM and in Run 2 
of Zhuge et al. (1994). Nevertheless, significant differences 
are visible: the maximum at t rj 2.5 ms in our model is 
10% higher and the width of the main peak about 40% 
narrower. Note that due to the difficulty of comparing the 
times for both models, the temporal positions of the first 
luminosity maxima were chosen to be the same and a nu- 
merical time shift of the maxima thus cannot be revealed. 
The largest difference of the gravitational wave luminosity 
appears at the secondary maximum at t ~ 4.6 ms after 

1 When the hydrodynamics code uses the relation of Eq. ^, 
heating means a change of the constant K in the polytropic 
equation of state, Eq. hi 



the actual merging event. In our simulations this maxi- 
mum is a factor of 2.5 higher than that shown in Zhuge et 
al. (1994). In addition, our Model ZCM displays a tertiary 
maximum at t » 7.3 ms, with the time intervals between 
the maxima being roughly 2.5 ms. The much more promi- 
nent secondary and tertiary maxima in our Model ZCM 
are due to the fact that the rotating merged object has a 
strongly quadrupolar structure and performs oscillations 
and pulsations: two revolving "subcores" can be seen in 
Fig. [j] which continue to emit gravitational waves for as 
long as they are present and rotate. 

The integrated energy emitted in gravitational waves 
is 2.0 • 10 52 erg which corresponds to 0.79% of Mc 2 , M 
being the mass of one neutron star. If one extracts the 
analogous values for the gravitational wave energy from 
Fig. 16b in Zhuge et al. (1994) for the time interval shown 
in Fig. |^, one obtains 0.8% of Mc 2 , in very good agreement 
with our Model ZCM, because the smaller peaks of the 
gravitational wave luminosity in the model of Zhuge et 
al. (1994) are compensated by their larger widths. 

For completeness we also show in Fig. |^ the gravita- 
tional wave forms for Model ZCM, which can be compared 
with Fig. 15 of Zhuge et al. (1994). After the maximum 
amplitude is reached, the wave form of model ZCM is 
damped by roughly a factor of two within 10 further pe- 
riods, while the damping factor is about 20 in the simu- 
lations of Zhuge et al. (1994). This difference reflects the 
motion of the two persisting "subcores" in our model. 

3.4- Gravitational wave spectrum 

The cumulative emission of gravitational wave energy as a 
function of frequency is shown in Fig. ||. In the upper part 
of each panel, the downward sloping straight line repre- 
sents the energy loss per unit frequency interval of a point- 
mass binary. The frequencies to the left of the vertical line 
correspond to the wave frequencies that are emitted before 
the start of the numerical modelling, i.e. for times t < 0. 
We produce a combined wave by using the quadrupole 
moments for a point-mass binary for t < and taking the 
numerically obtained quadrupole moments at t > (in 
analogy to Zhuge et al. 1994). The combined wave is then 
Fourier analysed up to the times given in the lower left 
corners of the panels of Fig. |[ 

At low frequencies the energy spectrum calculated for 
the combined wave fits very well to what is expected from 
the point-mass approximation. Although the merged ob- 
ject in the simulation radiates gravitational waves for a 
longer time than the point-masses, the very much smaller 
amplitudes associated with the coalescence of extended 
neutron stars result in less energy emitted at almost all 
higher frequencies. By Fourier transforming the signal of 
the combined wave until different times t > 0, we are able 
to roughly locate the moments when the peaks are pro- 
duced. 
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Fig. 5. Energy spectrum of gravitational waves emitted in Model ZCM. The times until which the Fourier transforms are 
performed are indicated in the lower left corner of the panels. They correspond to the snapshots of the density distributions 
shown in Fig. |l]. The straight, downward sloping line is the spectrum of a point-mass binary. The vertical lines indicate the 
frequency corresponding to the orbital frequency of two point masses at the initial distance of the neutron stars in our numerical 
model 



At around t w 2.3 ms the spectrum is structure- 
less. A peak at roughly / w 1.14 KHz starts to form 
at t « 3.24 ms and has saturated by t w 5.79 ms and 
shifted to around / w 1.22 KHz. A second peak around 
/ « 1.67 KHz is visible at this time but continues to grow 
until t « 8.7 ms and to shift to / « 1.79 KHz. The timing 
of the two peaks corresponds very well to the three max- 
ima of the gravitational wave luminosity (Fig. ||): the peak 
at / « 1.22 KHz is formed by the first outburst of grav- 
itational wave emission and the / rs 1.79 KHz peak by 
the second and third phases of strong emission. Fig. 17 in 



Zhugc et al. (1994) displays the corresponding wave spec- 
trum for their model. One notices that (a) we confirm 
their dip below 1 KHz, (b) we also see the slight rise above 
1 KHz, (c) they quote a value for a peak at 1.75 KHz which 
we indeed find at 1.79 KHz, and (d) their 1.75 KHz peak 
is significantly smaller than the one of our model. Thus 
the overall shape of the spectrum of our Model ZCM is 
similar to the spectrum of Run 2 of Zhuge et al. (1994). 
The visible differences, in particular the high frequency 
spectrum and the spectral maximum at / ks 1.75 KHz 
(point (d) above), are associated with the fact that due to 
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the quadrupolar deformation of the merged object in our 
simulation the emission of gravitational waves continues 
for a longer time than in Run 2 of Zhuge et al. (1994). This 
emission of waves at higher frequencies is an indication 
that the merged object becomes more compact with time. 
The gravitational waves emitted between 4 ms and 8 ms 
(Fig. ||) have a frequency of roughly 1.8 KHz, while those 
originating between 2 ms and 4 ms have around 1.2 KHz. 

We can only speculate why our Model ZCM wobbles 
around and oscillates for a longer time than the merger 
in the simulation of Zhuge et al. (1994). To calculate 
their model they used only 1024 SPH particles, which is 
at least an order of magnitude less than the number of 
zones resolving the central part of the merged object in 
our Model ZCM. Thus it seems possible that in Zhuge 
et al. (1994) a larger numerical viscosity of the SPH code 
damps out the quadrupolar pulsations more quickly than 
in our calculations. Another possible cause for the differ- 
ent oscillation behaviour of our Model ZCM compared to 
Zhuge et al. (1994) might be the differences of the phy- 
sical description, namely the lack of tidal deformations in 
the initial state and the implementation of backreaction 
terms in Model ZCM. The fact that Model ZCM does not 
initially have any tidal bulges and thus is not in equi- 
librium with the gravitational field will lead directly to 
oscillations of the neutron stars reflected e.g. in the values 
of the maximum density. This has been observed clearly in 
previous models, e.g. see Fig. 10 in Ruffert et al. (1996), 
in which one sees that the oscillations have a period of one 
dynamical time (i.e. the fundamental mode) and damp out 
in roughly 10 periods. Although the wobbling noticed in 
model ZCM might conceivably be triggered by this initial 
disequilibrium, the former is surely not sustained by the 
latter: the wobbles continue for a time over which the ini- 
tial oscillations are damped out. The effect of the different 
implementation of the gravitational wave backreaction is 
more difficult to assess, but in general it has a damping 
effect. Initially, when the neutron stars are well separated, 
the backreaction terms are responsible for the decay of the 
orbit, i.e. the spiral-in of the neutron stars. This phase is 
simulated with similar effect by Model ZCM and Zhuge et 
al. (1994), since without gravitational waves the neutron 
stars would circle each other practically indefinitely. How- 
ever, during the final plunge and the subsequent evolution 
Zhuge et al. (1994) do not include any backreaction. Once 
the neutron stars have come close enough, hydrodynamic 
effects alone dominate the merging process (cf. Rasio & 
Shapiro 1994 and references therein). Our Model ZCM is 
thus subjected to more damping during the final stages of 
the merging process than in Zhuge et al. (1994), so this 
again cannot explain the longer-term wobbling. 



4. Results of Model SNO 



4.1. Initial conditions 

In Model SNO we use the same mass and radius for the 
neutron stars and the same adiabatic exponent (r = 2) as 
Shibata et al. (1992) for their Model III. Thus, we place 
the two neutron stars with a mass of M = 1.4 Mq each 
and a radius of 9 km at an initial distance of 27 km. The 
two main differences between our run and the one of Shi- 
bata et al. (1992) are (1) the hydrodynamical integrator, 
and (2) the initial rotational state of the neutron stars 
and the correspondingly constructed equilibrium config- 
uration. Although the algorithm to calculate the grav- 
itational wave backreaction terms is the same in both 
works, we integrate the hydrodynamic equations with the 
PPM scheme, in contrast to Shibata et al. (1992) who cm- 
ployed a finite difference method that is not based on a 
Riemann-solver. In a number of previous models Naka- 
mura & Oohara (1991) considered neutron stars at rest in 
a corotating system which is equal to a solid-body type 
rotation in an inertial frame. Contrary to that, Model III 
of Shibata et al. (1992) was constructed as a "spinning" 
model, because the neutron stars were given spins in the 
corotating frame. Corresponding to this neutron star rota- 
tion Shibata et al. (1992) assumed axisymmetric neutron 
stars in rotational equilibrium initially, but did not solve 
for the equilibrium configuration in the common gravita- 
tional potential, arguing that tidal forces are much smaller 
than self-gravity. However, with the velocity distribution 
chosen by Shibata et al. (1992) the initial velocities of all 
parts of the neutron stars are collinear in an inertial frame, 
say parallel to the y-axis, and the neutron stars do not ro- 
tate. We therefore construct our neutron star models with 
the same velocity distrubution as unperturbed spherical 
poly tropes, whereby we also assume that tidal forces are 
small compared to self-gravity. We consider it to be more 
plausible that, if no spin is present in the inertial system, 
spin rotational deformation docs not need to be taken into 
account either. 

In their publication Shibata et al. (1992) show the den- 
sity distribution and flow pattern as functions of time for 
Model III as well as the gravitational wave forms and the 
gravitational wave luminosity. We shall therefore concen- 
trate on these quantities in the comparison of Model III 
of Shibata et al. (1992) and our Model SNO. 

4-2. Dynamical evolution 

Snapshots of the density distribution in the orbital plane 
together with the velocity field are shown at six instants 
in Fig. ||. The panels could, in principle, be compared 
with the panels of Fig. 3 of Shibata et al. (1992), respec- 
tively. In practice, however, the direct comparison of these 
contour plots is difficult, unfortunately, since Shibata et 
al. (1992) do not use contour values that are constant in 
time and, moreover, they plot arrows to indicate the flow 
field, but do not specify the actual magnitude of the fluid 
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Fig. 6. Cuts in the orbital plane of Model SNO at six instants in time showing the density contours together with the velocity 
field. The density contours are linearly spaced in intervals of 0.1p max starting at 0.01p m ax and are labeled in units of 10 14 g/cm 3 . 
The legend at the top right corner of each panel gives the scale of the velocity vectors and the time elapsed since the beginning 
of the simulation 



10 M. Ruffert et al.: Coalescing neutron stars - gravitational waves from polytropic models 

grav. wave luminosity, model SNO 




0.0 0.5 1.0 1.5 2.0 

time [ms] 



Fig. 9. The gravitational wave luminosity as a function of time 
for Model SNO (bold). The crosses represent the values taken 
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Fig. 8. Kinetic energy, internal energy, gravitational potential 
energy, and emitted gravitational wave energy as functions of 
time for Model SNO. The total energy contains all individual 
energies 



velocities. We used their prescription for the values of the 
density contours in Fig. [| 

One can recognize that initially, at t k, 0.4 ms and 
at t w 0.8 ms, our Model SNO is still in agreement with 
the one of Shibata et al. (1992). However, at t sa 1.2 ms 
two distinct minima appear between the cores of the neu- 
tron stars in Model SNO, while in Model III of Shibata 
et al. (1992) the most prominent minimum is at the cen- 
ter in between the neutron stars. These two minima of 
our Model SNO merge to form the central minimum at 
t Rj 1.6 ms, while two secondary minima are addition- 
ally present to the sides of the central minimum. Thus 
Model SNO at t « 1.6 ms resembles Model III, albeit at 
a different time, t & 1.2 ms. This indicates a shift in time 
by roughly At w 0.3 ms of the temporal evolution of both 
models that will become more clear in the next section. 
The three minima between the neutron star cores merge 
again, forming one elongated string by t w 1.8 ms. The 
structure of the merged object with two sub-cores remains 
present until the end of our simulations at t ~ 2 ms, while 
the density valley between becomes less prominent. 

A quantity better suited for comparison is the max- 
imum density as a function of time, the values of which 
can be found inserted in the panels of Fig. 3 of Shibata 
et al. (1992). In Fig. [7] these values are displayed at eight 
points in time, together with the maximum density for 
our Model SNO. The maximum density of Model III by 
Shibata et al. (1992) peaks at around t s=s 0.4 ms which 
is different from our Model SNO. It is hard to explain 
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Fig. 10. The gravitational wave forms, h+ and h x , for 
Model SNO as observed on the z-axis at 10 Mpc distance 

this difference of the density evolution without having 
more detailed information about Model III of Shibata et 
al. (1992). At the time the density maximum occurs in the 
latter model, the neutron stars are still clearly separated, 
in both their Model III (Fig. 3c in Shibata et al. 1992) 
and in our Model SNO (Fig. g). There is an indication 
that at t « 0.4 ms the density distribution in the neu- 
tron stars of Shibata et al. (1992) develops a "crack" or 
"discontinuity". The "crack" seems to emanate from the 
center and propagate along the direction of the x- and y- 
axes in both neutron stars. We speculate that this might 
cause the higher density maximum at t = 0.4 ms in the 
model of Shibata et al. (1992). The reason for the differ- 
ences in the evolution of the maximum density and the 
neutron star separation seems not to be associated with 
the initial shape and structure of the neutron stars due 
to the different assumptions about the neutron star spins. 
Rather than differences of the mass distributions, the dif- 
ferent momenta of the stars might be responsible for the 
discrepant evolution. 

In Fig. H we show the temporal evolution of some in- 
tegral energy quantities: the potential, kinetic, internal 
energies as well as the energy in gravitational waves and 
the total sum of all energies. The total energy is conserved 
to 1% compared to the maximum potential energy. The 
more compact structure of the merged object is reflected 
in the deeper potential at the end of the simulation. Note 
the slight oscillation present especially in the potential 
and kinetic energies which shows the wobbling of the sub- 
cores. The internal energy hardly changes indicating that 
the merging process proceeds rather softly without strong 
shocks. 



4-3. Gravitational wave forms and luminosity 

Figure |9| displays the gravitational wave luminosities as 
functions of time for Model SNO and for Model III of Shi- 
bata et al. (1992; Fig. 6b). The maximum of the emission 
(which is located in the model of Shibata et al. (1992) at 
t ss 0.9 ms) is shifted in time by At « 0.3 ms, but agrees 
in height to within 10%. Note that model SNO follows 
the point-mass binary approximation for roughly 0.2 ms 
longer than Model III. We assume that (a) this initial de- 
viation of Model III is primarily responsible for the sub- 
sequent time shift visible in the density contours and the 
gravitational wave luminosity, and that (b) the initial de- 
viation is due to the differing initial rotational states of 
the neutron stars. Since the separation at the beginning 
of the simulations is very close to the dynamic stability 
limit (cf. Lai et al., 1994, and references cited therein) a 
small difference in density distribution can cause a sub- 
stantial time lag: if the neutron stars are initially further 
separated than the stability limit, gravitational wave emis- 
sion has to continue for some time to drive the decay of 
the orbit before the dynamical instability can set in. 

For completeness we also show in Fig. [to] the gravita- 
tional wave forms for Model SNO which can be compared 
with Fig. 6a of Shibata et al. (1992). Due to the time shift 
mentioned above, the signal of our Model SNO is roughly 
one wave period shorter at the end of the simulation than 
Model III. 



5. Results of Model RJS 

5.1. Initial conditions 

Model RJS is calculated with neutron stars which have 
the same mass, radius, and central density as those in 
Model A64 of Ruffert et al. (1996). Using the Lanc-Emden 
equation and the mass-radius relation, this was achieved 
by appropriately adjusting the polytropic constant K and 
the adiabatic exponent T (cf. Table |l|) for the equation 
of state. Placing the neutron stars at the same initial dis- 
tance of 42 km and giving them the same velocity distribu- 
tion as in as Model A64, we are left with a model differing 
only in the employed equation of state: Model RJS uses 
a polytopic equation of state according to Eq. ||, while 
Model A64 was calculated with the equation of state of 
Lattimer & Swesty (1991) and included the effects due to 
neutrino emission. 

Out of the variety of quantities shown and discussed 
in Ruffert et al. (1996), we decided to compare the evolu- 
tion of the mass distribution, the separation of the density 
maxima as a function of time to show possible differences 
in the dynamics, and the luminosity, wave form, and en- 
ergy spectrum of the gravitational wave emission. 
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Fig. 11. Cuts in the orbital plane of Model RJS at four instants in time showing the density contours together with the velocity 
field. The density contours are logarithmically spaced with intervals of 0.5 dex. The density is measured in units of g/cm 3 . The 
bold contours are labeled with their respective values. The legend at the top right corner of each panel gives the scale of the 
velocity vectors and the time elapsed since the beginning of the simulation 



5.2. Dynamical evolution 

The contour plots in Fig. [ll] show the dynamical evolu- 
tion. The mass distribution is visualized by density con- 
tours with arrows superimposed to indicate the flow ve- 
locities. These plots can directly be compared with the 
Figs. 4c, 4e, 5a, and 5c of Ruffert et al. (1996), respec- 
tively. Initially, until about t w 1.7 ms, only minor dif- 
ferences occur. At t w 1.7 ms the neutron star surfaces 
towards the downstream sides (at the positions i«0 km, 
y w ±20 km) are more extended in Model RJS (Fig. [Tl|b) 
than in Model A64 of Ruffert et al. (1996) (see Fig. 4e 
there). This effect also can be seen in Fig. pd| : where the 



contour for the density value of 10 11 g/cm 3 extends fur- 
ther out than the marginal "spiral arm" visible in Fig. 5a 
in Ruffert et al. (1996). 

Compared to the models computed with the equa- 
tion of state of Lattimer & Swesty (1991) in Ruffert et 
al. (1996), all three models discussed here, ZCM, SNO, 
and RJS, have a less steep density decline at the neutron 
star surface and develop a more extended disk after merg- 
ing. This general phenomenon is explained by the stiffness 
of the polytropic equations of state which have adiabatic 
exponents between 2 and 2.32. This is much stiffer than 
the equation of state of Lattimer & Swesty (1991) in the 
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Fig. 12. The separation of the density maxima of the two 
neutron stars as a function of time for Model RJS (bold) and 
for Model A64 (dotted) of Ruffert et al. (1996) 



Fig. 13. The maximum density on the grid as a function of 
time for Model RJS (bold) and for Model A64 (dotted) of Ruf- 
fert et al. (1996) 



subnuclear regime, where T ss 1.33 to 1.66. An increase 
of the internal energy, e.g. by the dissipation of kinetic 
energy in shocks or via shear effects, which is important 
in the regions close to the surface, therefore, leads to a 
strong pressure increase and thus to an inflation of the 
stellar layers involved. 

In the core region of the merged object, Figs. |lT| c 
and |TT[ d reveal more fine structure than Model A64 of 
Ruffert et al. (1996; Fig. 5a and 5c) which indicates a 
stronger damping of high-frequency oscillation modes in 
the latter model, caused by the bulk viscosity of the phys- 
ical equation of state used there. At time t w 4.5 ms 
both simulations show a rather similar overall structure 
of the merger (compare Fig. [Tl| d to Fig. 5c in Ruffert et 
al. 1996). In both cases the innermost density contour 
(p » 3 ■ 10 14 g/cm 3 ) is dumb-bell shaped and the underly- 
ing double-core structure means that the merging process 
is still going on. Also during this phase the exterior layers 
at subnuclear densities are significantly more extended in 
Model RJS due to the stiffer polytropic equation of state. 

On the whole, the coalescence in Model RJS pro- 
ceeds very similarly to Model A64 of Ruffert et al. (1996). 
Figs. [l2| and |l3| show the separation of the density max- 
ima and the value of the maximum density as functions 
of time for the two models. In both models the neutron 
star density maxima merge within 2 ms and then sepa- 
rate again for a short while when the double-core struc- 
ture forms (between 2.5 ms and about 3.5 to 4 ms) be- 
cause of the large angular momentum of the merged ob- 
ject. Model RJS comes to "rest" faster (by roughly 0.5 ms) 
than Model A64 which has separated density maxima 



from 2.5 ms to 4.5 ms. This difference might be caused by 
a larger angular momentum transport outwards through 
the more extended disk of Model RJS. This hypothesis is 
strongly supported by the outwardly directed velocity field 
which can be seen in Fig. [TT|d, indicating mass and angular 
momentum loss accross the grid boundaries. Notice also 
that due to the differences of the subnuclear equation of 
state the transition from the central core of the merged 
object to the disk region is less sharp in Model RJS (com- 
pare Fig. [ll|b with Fig. 5c in Ruffert et al. 1996). 

Before the neutron stars merge, their central density 
decreases slightly by about 5% (see Fig. [l3]) and during 
the coalescence the maximum density increases by about 
20% of the initial value. The overall dynamical behaviour, 
reflected by the separation of the density maxima and the 
value of the maximum density, reveals much less difference 
between the models than can be seen from the density 
contours. The similar dynamics can easily be understood 
because the construction of the initial neutron stars and 
the adjustment of the parameters K and T in the poly- 
tropic equation of state Eq. |l| should ensure very similar 
behaviour (e.g response to compression) of the polytropic 
equation of state and the Lattimer & Swesty (1991) equa- 
tion of state in the nuclear and supranuclear regimes. The 
differences occur in the outer regions of the neutron stars 
at densities below the nuclear matter density where the 
polytropic equation of state for the chosen value of T is 
much stiffer than the equation of state of Lattimer & 
Swesty (1991). This difference mainly concerns the evo- 
lution of the outer layers and thus the formation, struc- 
ture and properties of the disk that forms during merging. 
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Fig. 14. The gravitational wave luminosity as a function of 
time for Model RJS (bold) and Model A64 (dotted) of Ruffert 
et al. (1996) 



Since, however, only a comparatively small fraction of the 
neutron star mass has densities less than nuclear matter 
density, the overall dynamical evolution of the bulk of the 
matter is unaffected by the difference of the low-density 
equation of state. From the gradual increase of the density 
maximum we conclude that the merging is as smooth as 
in Model ZCM and not so violent as in Model SNO where 
the density increases rapidly during the final coalescence 
(see Fig. (jj). 

5. 3. Gravitational wave forms and luminosity 

The gravitational wave luminosity is a very sensitive in- 
dicator of differences in the overall dynamical evolution 
and the mass distributions, since it involves the third 
time derivative of the quadrupole moment. One recognizes 
from Fig. 14 that the main maximum of the gravitational 
wave luminosity is lower by 20% and slightly narrower in 
Model RJS compared to Model A64. Roughly the same 
statements apply to the second maximum. These differ- 
ences between the two models are concordant with the 
temporal behaviour of the separation of the density max- 
ima in both models (see Fig. |l2]), but are much smaller 
than the differences between Models A64, B64, and C64 
of Ruffert et al. (1996). 

The wave forms (Fig. |l^) show only minor differences: 
both the amplitudes and the phases of the signals coincide 
well. Slight discrepancies develop after t « 2 ms which cor- 
responds to the time after the first merging of the density 
maxima has occurred (Fig. 12) and the gravitational wave 
emission reaches a maximum (Fig. |l4| ). Once the separate 
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Fig. 15. The gravitational wave forms, h+ and h x , for 
Model RJS (bold) and Model A64 (dotted) of Ruffert et 
al. (1996) 



density maxima have vanished, the merged configuration 
has loses its strongly quadrupolar deformation and there- 
fore radiates gravitational waves less efficiently. 

In agreement with the temporal displacement of the 
evolution visible in Figs. [l2], [r], and [l4|, a small phase 
shift of the gravitational wave forms after t i=s 2 ms is 
caused by the slightly faster rotation of Model A64. The 
higher angular momentum of the merged configuration 
(see above) together with the more compact structure of 
the layers at intermediate densities (between 10 12 g/cm 3 
and 10 14 g/cm 3 ) due to the less stiff equation of state, is 
also the reason for somewhat larger gravitational wave 
amplitudes and the enhanced gravitational wave lumi- 
nosity in case of Model A64 at times later than about 
t « 3 ms. The total energy emitted in gravitational waves 
in Model RJS is 2.8 • 10 52 erg, roughly 20% less than the 
3.4 • 10 52 erg emitted by Model A64 within 5 ms (taken 
from Fig. 23 of Ruffert et al. 1996). 

5.4. Gravitational wave spectrum 

The close resemblance of the gravitational wave forms 
(Fig. [l5|) is mirrored in the gravitational wave spec- 
tra which share the main features in the two mod- 
els RJS and A64 from Ruffert et al. (1996). The spectra 
of Model RJS can be found in Fig. |l6| and should be com- 
pared with Fig. 28 in Ruffert et al. (1996). They contain 
the same kind of information as already described in the 
context of Fig. |[ i.e. the frequency distribution of the 
gravitational wave energy emitted until the times shown 
in the lower left corners of the plots, together with a line 
representing the energy loss per unit frequency interval of 
a point-mass binary. 
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At time t w 3.2 ms the spectra look very similar and 
Models RJS and A64 both have a spectral minimum at 
/ ~ 1 KHz and spectral maxima at / w 1.3 KHz and / w 
2.0 KHz. Even for frequencies / > 2.0 KHz the decline 
of the spectra is pretty much the same. The major visi- 
ble differences are the strengths of the second and of the 
third spectral maximum at / « 2.4 KHz. While the max- 
imum at / w 2 KHz is more pronounced in Model RJS, 
the peak at / « 2.4 KHz is higher than in Model A64. In 
both models the / « 1.3 KHz maximum is mainly gener- 
ated before t ~ 3.2 ms, while the / « 2 KHz maximum 
grows strongly after t « 3.2 ms. The / ~ 2 KHz maximum 
tends to shift to f ks 1.8 KHz (more clearly in Model RJS) 
between t w 3.2 ms and t w 4.64 ms and becomes slightly 
more prominent in Model RJS than in Model A64. Af- 
ter t f=s 4.6 ms the spectral peaks at frequencies between 
about 2 KHz and 2.4 KHz have merged into one struc- 
ture in Model A64, while in Model RJS there are still two 
sharp, well separated maxima and the / w 2.4 KHz fea- 
ture has not gained much strength any more. We think 
that these results can again be explained by the stronger 
pulsational and vibrational activity of Model A64 after 
merging which adds to the power of the gravitational wave 
emission at frequencies around and above / 2 KHz. 

6. Discussion and summary 

The hydrodynamic simulations presented in this work fol- 
low the dynamical evolution and the gravitational wave 
emission of two merging neutron stars modeled as poly- 
tropes. We compared our results with three models pub- 
lished earlier and tried to obtain information about uncer- 
tainties and differences caused by the numerical scheme 
employed for integrating the hydrodynamic equations, 
by the numerical resolution and setup of the initial 
model, and by the use of a polytropic equation of state 
instead of the physical equation of state of Lattimcr 
& Swesty (1991). We chose the following three models 
for comparison: (a) "Run 2" of Zhuge et al. (1994), (b) 
"Model III" of Shibataet al. (1992), and (c) "Model A64" 
of Ruffert et al. (1996). We attempted to construct the ini- 
tial configurations as similar to these models as possible 
and also used polytropic equations of state which were 
the same or which guaranteed to reproduce the properties 
(mass, central density, radius) of the neutron stars. 

6.1. Model ZCM 

Run 2 of Zhuge et al. (1994) is different from our 
Model ZCM in (a) one numerical aspect: Zhuge et 
al. (1994) used an SPH algorithm while we employed the 
grid-based PPM scheme and (b) two physical aspects: (i) 
the gravitational wave backreaction is prescribed in an ad 
hoc way by Zhuge et al. (1994), because it is included in 
the computation when the neutron stars are still separate 
but switched off during the merging process; (ii) both com- 



putations start with different initial separations and tidal 
deformations. Our Model ZCM seems to retain a "double 
core" structure and thus a large quadrupolar deformation 
during the merging for a longer time than Run 2. This 
has two consequences: On the one hand, the secondary 
and tertiary maxima of the gravitational wave luminos- 
ity are significantly higher in our calculation, and on the 
other hand, the peak in the gravitational wave spectrum 
at a frequency of / ss 1.8 KHz is more prominent, too 
(a first maximum is present at about 1.3 KHz). Both ef- 
fects are caused by a more extended post-coalescence pe- 
riod in Model ZCM where the merger performs "ringing" 
and "wobbling" motions and emits more power in grav- 
itational waves at frequencies around / « 1.8 KHz. We 
find a local minimum of the gravitational wave spectrum 
at a frequency of about / ps 1.5 KHz in good agreement 
with Run 2 of Zhuge et al. (1994), although this mini- 
mum is not quite as deep in our calculation. The inte- 
grated gravitational wave energy in both models is also 
nearly the same, implying that the stronger emission of 
our Model ZCM at high frequencies is accompanied by 
a deficit at lower frequencies, which we indeed find in 
the band / ps 0.7-1.1 KHz. We suppose that these dif- 
ferences are mainly caused by the larger numerical viscos- 
ity of the SPH code and the comparatively poor resolu- 
tion with only 1024 particles in the computation of Zhuge 
et al. (1994). Although the two main physical differences 
(initial tidal deformations, backreaction implementation) 
could in principle be also partly responsible for the wob- 
bling, we assume their effect to be small: (i) the initial dise- 
quilibrium due to tidal deformations has been shown (Ruf- 
fert et al. 1996) to be damped out within a time shorter 
than the ringing, and (ii) the gravitational backreaction 
has a damping effect in the longer-term merging process. 

6.2. Model SNO 

Shibata et al. (1992) performed their calculations with a 
finite difference grid-based explicit Eulerian code just as 
we do, albeit with a different hydrodynamical integrator, 
a second order advection using a method proposed by 
LcBlanc (sec Appendix A of Oohara & Nakamura 1989), 
which does not make use of a Riemann solver. They in- 
corporated the gravitational wave backreaction in a way 
conceptually equivalent to our treatment. Their grid had 
121 3 zones. To solve the Poisson equation, the Incomplete 
Cholcsky decomposition and Conjugate Gradient (ICCG) 
method by Meijerink & van der Vorst (1977) was used. 
Contrary to this we implemented a fast Fourier method 
and our grids had a size of 128x128 (and 32 zones verti- 
cal to the orbital plane). However, long before any hydro- 
dynamical interaction between the neutron stars occur, 
at t ps 0.4 ms when the neutron stars just barely touch 
each other, the maximum density on the grid showed an 
extreme value in the calculations of Shibata et al. (1992) 
which we were unable to reproduce in our Model SNO. 
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Fig. 16. Energy spectrum of gravitational waves emitted in 
Model RJS. The times until which the Fourier transforms are 
performed are indicated in the lower left corners of the panels. 
They correspond to the snapshots of the density distributions 
shown in Fig. 28 of Ruffert et al. (1996). The straight, down- 
ward sloping line is the spectrum of a point-mass binary. The 
vertical lines indicate the frequency corresponding to the or- 
bital frequency of two point masses at the initial distance of 
the neutron stars in our numerical model 
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This difference is difficult to explain and we can only spec- 
ulate that it might be associated with the choice of Shibata 
et al. (1992) to use configurations in rotational equilib- 
rium for the two neutron stars which were given spins in 
the frame that corotates with their orbital revolution. Be- 
cause these spins of the stars were absent for an observer 
in an inertial frame and because Shibata et al. (1992) did 
not construct equilibrium configurations in the common 
gravitational field of the two stars (which would make 
more sense), we did not follow their choice but started 
our simulations as usual with spherical neutron star mod- 
els. This obvious discrepancy causes an initial delay of the 
merging by At « 0.3 ms in our model SNO compared to 
Model III of Shibata et al. (1992). Nevertheless, the grav- 
itational wave luminosity amplitude of their model and of 
our simulation coincide surprisingly well to within 10%. 



6.3. Model RJS 

In a third comparative study we performed a simulation, 
Model RJS, with initial conditions (neutron star masses, 
radii, central densities, initial separation, and velocities) 
like those used in Model A64 of Ruffert ct al. (1996). The 
only major difference was that instead of the physcial 
equation of state of Lattimer & Swesty (1991) — which 
also allowed us to include neutrino physics in the com- 
putation — we used a polytropic equation of state in 
Model RJS with the initial value of the parameter K and 
the constant adiabatic index T adjusted such that the ini- 
tial neutron star properties are reproduced. 

We found the overall dynamical evolution during the 
coalescence, e.g. timescale of merging, formation of spiral- 
arm-like structures of matter spun off the neutron star 
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surfaces by tidal and centrifugal forces, post-merging os- 
cillations and pulsations, to be very similar in both mod- 
els. This is particularly clearly seen in the behaviour of the 
neutron star separation and of the maximum density on 
the grid as functions of time. Correspondingly, the grav- 
itational wave luminosity, gravitational wave form, and 
gravitational wave spectrum share the main structural fea- 
tures. The peak gravitational wave luminosity and the to- 
tal energy emitted in gravitational waves are only about 
20% smaller in the polytropic model. The gravitational 
wave forms are very similar but have the tendency to de- 
velop a slight phase shift after the two neutron stars have 
merged into one object. These minor differences reflect the 
fact that the properties of the Lattimer & Swesty (1991) 
equation of state in the nuclear and supranuclear regimes 
- which are most important for the neutron star struc- 
ture as well as for the gravitational wave signature dur- 
ing the merger — are sufficiently well reproduced. This is 
ensured by our determination of the equation of state pa- 
rameters K and L by prescribing the same mass, central 
density and radii of the neutron star models constructed 
with the polytropic equation of state and with the Lat- 
timer & Swesty (1991) equation of state. 

The minor differences of both calculations, the roughly 
20% smaller peak emission of gravitational waves in the 
polytropic case and the small phase shift of the grav- 
itational wave forms, are mainly caused by differences 
of the equations of state in the subnuclear regime (p 
10 g/cm 3 ). In this regime in particular, the description 
of the properties of neutron star matter with an equa- 
tion of state with a uniformly chosen adiabatic exponent T 
is inadequate. With the employed value of T = 2.319 
the equation of state at subnuclear densities is by far 
too stiff. As a result, the neutron star surface layers in 
Model RJS were more extended and the density decline 
less sharp compared to Model A64 in Ruffert et al. (1996). 
Morerover, during merging the neutron star surfaces came 
into contact earlier in Model RJS which moderated the 
dynamical interaction of the stars and reduced the lumi- 
nosity outburst of gravitational waves during this phase. 
After merging, the overestimated stiffness of the equation 
of state led to a slightly less compact central body of the 
merger in Model RJS and to a significantly more extended 
disk structure surrounding the central, massive body than 
in Model A64. In particular, the properties of this disk, its 
mass, radial extent, and rotational state, therefore require 
the use of a proper physical equation of state during the 
simulations of neutron star merging. 
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